#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup2.R")

#讀檔案
dta <- read.csv(file = "../Data/HEXACO_HS.csv", 
                  header = TRUE, stringsAsFactors = TRUE)

#檢視資料結構
#程式報表2.1
str(dta)

#設定類別變項的順序
dta <- dta |> 
  dplyr::mutate(母親教育程度 = forcats::fct_relevel(母親教育程度, 
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")),
                父親教育程度 = forcats::fct_relevel(父親教育程度, 
                                              c("小學或不識字", 
                                                "國中", 
                                                "高中", 
                                                "大學或專科", 
                                                "研究所以上")))

#瞭解 ggplot 啟動後的設定規劃
#程式報表2.2
ggplot() |> attributes()

#ggplot 是一個一個圖層（layer）疊上去
#底下以散佈圖為例，呈現每一步驟的結果
#第一步驟設定圖的框架，注意圖的X軸與Y軸
#圖2.1
g0 <- ggplot(data = dta, 
             aes(x=社會退縮, y=攻擊行為)) 
g0

#在前一步驟結果（圖的框架）上加入點
#圖2.2
g1 <- g0 + geom_point(alpha=.2) 
g1

#在前一步驟結果上加入橢圓（資料的95%區間）與局部迴歸線
#圖2.3
g2 <- g1 + stat_smooth(method='lm', 
                       formula = y ~ x,
                       se=FALSE, 
                       linewidth = 0.5) 
g2

#在前一步驟結果上要求以變項（性別）區分顏色
#圖2.4
g3 <- g2 + aes(color=性別) 
g3

#在前一步驟結果上設定 x 軸與 y 軸刻度
#圖2.5
g4 <- g3 + 
  scale_y_continuous(breaks=seq(0, 25, by=5)) +
  scale_x_continuous(breaks=seq(0, 12, by=2))
g4

#在前一步驟結果上要求以變項（父母教育）分面（facet）
#圖2.6
g5 <- g4 + facet_wrap(vars(母親教育程度),nrow=1)
g5

#在前一步驟結果上加上X軸、Y軸的標示，以及整個圖形的標題
#圖2.7
g6 <- g5 + labs(x='社會退縮分數',
                y='攻擊行為分數',
                title='散布圖：攻擊行為與社會退縮')
g6

#在前一步驟結果上改變主題，並要求圖示位置
#圖2.8
g7 <- g6 + theme_minimal() + 
  theme(legend.position='top')
g7

#前面步驟可以一次執行
ggplot(data = dta, aes(x=社會退縮, y=攻擊行為)) +
   geom_point(alpha=.2) +
 stat_smooth(method='lm', formula = y ~ x,
       se=FALSE, linewidth = 0.5)+
 aes(color=性別) +
 scale_y_continuous(breaks=seq(0, 25, by=5)) +
 scale_x_continuous(breaks=seq(0, 12, by=2)) +
 facet_wrap(vars(母親教育程度),nrow=1)+
 labs(x='社會退縮分數', y='攻擊行為分數',
     title='散布圖：攻擊行為與社會退縮') +
 theme_minimal() + 
 theme(legend.position='top')

#設定後面的 ggplot 繪圖以 theme_minimal 為預設值
ggplot2::theme_set(theme_minimal())

#以 ggplot 直接繪製統計摘要結果（連續資料平均數與標準誤） 
#圖2.9
ggplot(data=dta, 
     aes(x=母親教育程度, y=誠實.謙遜, color=性別)) +
  stat_summary(fun.data = "mean_cl_boot", 
               position=position_dodge(.2)) +
  stat_summary(aes(group=性別), fun = mean, 
               geom="line",
               position=position_dodge(.2))+
  scale_color_grey(end=.7)+
  labs(y='誠實.謙遜平均分數',
       x='母親教育程度',
       title='不同性別的誠實.謙遜平均跟母親教育程度的關係',
       caption="來源: 許功餘")+
  theme(legend.position='top')

#程式報表2.3
dta |> 
  dplyr::group_by(性別, 母親教育程度) |>
  dplyr::reframe( 誠謙平均 = mean(誠實.謙遜),
                  誠謙標準誤 = sd(誠實.謙遜)/sqrt(n()),
                  誠謙平均下界 = 誠謙平均 - 1.96*誠謙標準誤,
                  誠謙平均上界 = 誠謙平均 + 1.96*誠謙標準誤)

#以 ggplot 直接繪製統計摘要結果（類別資料的百分比） 
#圖2.10
ggplot(dta, 
       aes(x=母親教育程度, group=父親教育程度)) + 
  geom_bar(aes(y=after_stat(prop), 
          fill = factor(after_stat(x))),
          width = .2) + 
  scale_y_continuous(labels=scales::percent) +
  scale_fill_grey()+
  coord_flip()+
  labs(y="百分比",
       title="門當戶對:父母親教育程度") +
  facet_wrap(vars(父親教育程度), ncol=1)+
  theme(legend.position="none")

#程式報表2.4
dta |> 
  dplyr::select(母親教育程度, 父親教育程度) |>
  gtsummary::tbl_cross(percent=c("column"))

#以 ggplot 直接繪製統計摘要結果（連續資料平均數與標準誤） 
#圖2.11
ggplot(data=dta, 
     aes(x=母親教育程度, y=誠實.謙遜)) +
  stat_summary(fun.data = "mean_cl_boot") +
  scale_color_grey(end=.7)+
  facet_wrap(vars(父親教育程度), nrow=1)+
  labs(y='誠實.謙遜平均分數',
       x='母親教育程度',
       title='不同父親教育程度的誠實.謙遜平均跟母親教育程度的關係',
       caption="來源: 許功餘")+
  theme(legend.position='top',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

#多變量圖形 
#圖2.12,dpi=300
dta |> tidyr::pivot_longer(cols=4:9, 
                           names_to = '人格維度',
                           values_to = '分數') |>
  ggplot()+
  aes(x=違反規定, y=分數, color=性別)+
  geom_point(size=rel(.5))+
  stat_smooth(method='lm', 
              formula = y ~ x,
              se=F, linewidth=.5)+
  facet_grid(人格維度 ~ 母親教育程度) +
  scale_color_grey(start=.1, end=.6)+
  labs(x="違反規定分數",
       y="人格維度分數")+
  theme(legend.position='top')

#同時繪製兩個變項的散布圖，以及各自的邊際分布
#圖2.13
p <- ggplot(dta, 
            aes(x=誠實.謙遜, y=違反規定)) +
  geom_point(shape=21, alpha=.5) +
  stat_ellipse() +
  geom_vline(xintercept=mean(dta$誠實.謙遜), col="gray") +
  geom_hline(yintercept=mean(dta$違反規定), col="gray") +
  stat_smooth(method="lm", 
              formula = y ~ x,
              linewidth=.7,
              linetype="dotted",
              se=FALSE,
              alpha=.5,
              col=1) +
  labs(y="違反規定分數", 
       x="誠實.謙遜分數") 
pacman::p_load(ggExtra, KernSmooth)
ggMarginal(p, 
           type="histogram",
           xparams=list(binwidth=dpih(dta$誠實.謙遜),
                          fill="gray90"),
           yparams=list(binwidth=dpih(dta$違反規定),
                          fill="gray90"))

#以 tidy 方式整理後繪製
#圖2.14
lm(誠實.謙遜 ~ 性別 + 母親教育程度 + 性別:母親教育程度, data=dta) |>
  broom::tidy(conf.int=TRUE) |> 
  dplyr::slice(-1) |>
  ggplot() +
  aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0) +
  geom_errorbarh()+
  geom_point() +
  geom_vline(xintercept = 0, linetype='dotted', col='gray') 

#繪製模型分析所得參數
#以套件GGally整理後繪製
#圖2.15
lm(誠實.謙遜 ~ 性別 + 母親教育程度 + 性別:母親教育程度, data=dta) %>%
  GGally::ggcoef(.,  exclude_intercept=TRUE,
                 sort=NULL, 
                 na.rm=TRUE) +
  labs(x="估計值",
       y="變項")

#繪圖檢視模型與資料的配適性
#這邊先做性別與父母教育對數學成績的二因子 ANOVA
#程式報表2.5
lm(誠實.謙遜 ~ 性別*母親教育程度, data=dta) |> anova()

#繪製模型分析所得參數，並與實際資料對照
#先彙整資料
#程式報表2.6
ef_m1 <- lm(誠實.謙遜 ~ 性別 + 母親教育程度, data=dta) %>%
  ggeffects::ggpredict(., terms=c("母親教育程度", "性別")) |> 
  as.data.frame()

ef_m1

#繪製模型分析所得參數，並與實際資料對照
#注意 y 被設定四次
#圖2.16
ggplot() +
   stat_summary(data=dta, 
                aes(x=母親教育程度, y=誠實.謙遜, color=性別), 
                fun.data = "mean_cl_boot", 
                size=rel(.3),
                position=position_dodge(.2))+
   geom_line(data=ef_m1, aes(x=x, y=predicted, 
                             group=group, 
                             color=group))+
   geom_line(data=ef_m1, aes(x=x, y=conf.low, 
                             col=group, group=group),
                             linetype='dotted')+
   geom_line(data=ef_m1, aes(x=x, y=conf.high, 
                             col=group, group=group),
                             linetype='dotted') +
   scale_color_grey(end=.6)+
   labs(y='誠實.謙遜平均分數',
        x='母親教育程度',
        title='性別跟母親教育對誠實.謙遜平均的效果')+
   guides(color=guide_legend(title="性別"))+
   theme(legend.position='top')